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1.  FORWARD 

The  goal  of  this  project  has  been  to  develop  computational  tools  for  the  simulation  of 
candidate  airdrop  systems  for  the  New  World  Vistas  (NWV)  Precision  Airdrop  (PAD) 
program.  To  accomplish  this  goal,  researchers  from  the  University  of  Connecticut,  Rice 
University,  and  the  US  Army  Soldier  Systems  Center  at  Natick  have  collaborated  in 
developing  a  coupled  Computational  Fluid  Dynamics  (CFD)  and  Structural  Dynamics 
(CSD)  program  to  simulate  three-dimensional,  transient  Fluid-Structure  Interaction  (FSI) 
phenomena  for  airdrop  systems. 

The  work  performed  under  this  project  consists  of  (1)  development  of  new  CFD  methods, 
(2)  development  of  new  CSD  methods,  (3)  development  of  new  methods  for  coupled  FSI 
simulations,  (4)  verification  of  the  FSI  model,  and  (5)  simulation  of  candidate  NWV 
airdrop  systems.  The  FSI  simulations  require  large  scale,  nonlinear,  transient  finite 
element  models  for  the  parachute  system  and  surrounding  airflow  and  therefore  are 
computationally  intensive.  To  address  this  difficulty,  parallel  computational  techniques 
have  been  developed  for  the  FSI  model. 

The  CFD  model  is  based  on  the  Deforming-Spatial-Domain/Stabilized-Space-Time 
(DSD/SST)  formulation  developed  by  the  Team  for  Advanced  Flow  Simulation  and 
Modeling  (T*AFSM)  at  Rice  University.  The  CSD  model  is  based  on  the  Total 
Lagrange  formulation  for  geometrically  nonlinear  dynamic  behavior  of  flexible 
membranes  and  cables  developed  by  the  University  of  Connecticut.  The  coupled  FSI 
method  was  developed  by  Natick  and  Rice  University  in  collaboration  with  the 
University  of  Connecticut. 
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4.  PROBLEM  STATEMENT 


The  primary  goal  of  this  project  has  been  to  develop  computational  tools  for  the 
simulation  of  candidate  airdrop  systems  for  the  New  World  Vistas  (NWV)  Precision 
Airdrop  (PAD)  program.  The  computational  methods  require  the  coupling  of  a 
Computational  Structural  Dynamics  (CSD)  model  for  the  parachute  system  with  a 
Computational  Fluid  Dynamics  (CFD)  model  for  the  surrounding  airflow  to  perform 
Fluid-Structure  Interaction  (FSI)  simulations.  To  accomplish  this  goal,  the  following 
specific  problems  in  each  area  must  be  addressed: 

Computational  Structural  Dynamics 

•  A  robust  structural  model  is  required  to  simulate  the  transient,  geometrically 
nonlinear  behavior  of  thin  membrane  and  cable  structures.  Parachute  systems 
undergo  highly  transient  behavior  with  large  relative  motion.  Robust  algorithms 
are  required  to  maintain  stability  of  the  numerical  solution  and  accurately 
simulate  this  complex  behavior. 

•  A  rigorous  method  is  needed  to  model  the  loss  of  tension  in  fabric  membrane 
structures  such  as  parachute  canopies.  Fabric  structures  do  not  support 
compressive  stresses  and  undergo  “wrinkling”  at  the  onset  of  compression.  Since 
the  behavior  of  a  structure  will  differ  dramatically  depending  on  whether  there  is 
compression  or  no-compression,  it  is  necessary  to  account  for  wrinkling  in 
parachute  systems. 

•  Special  CSD  models  are  needed  to  model  Pneumatic  Muscle  Actuators  (PMAs) 
for  simulation  of  NWV  candidate  systems.  PMAs  are  the  primary  control  device 
being  used  in  the  NWV  PAD  program.  PMAs  utilize  a  braided  fiber  cylinder  that 
undergoes  large  relative  motions  under  rapid  pressurization.  A  geometrically 
nonlinear  formulation  for  PMAs  is  required. 

•  Stand-alone  CSD  simulations  are  needed  to  verify  the  structural  model  and  to 
provide  a  preliminary  modeling  capability  for  candidate  NWV  systems.  Stand¬ 
alone  CSD  simulations  with  approximate  fluid  loads  (prescribed  pressure  and 
velocity  dependent  drag) 

Computational  Fluid  Dynamics 

•  A  robust  fluid  dynamics  model  is  required  to  simulate  flow  fields  that  involve 
moving  boundaries  and  interface.  For  the  CFD  model,  it  is  necessary  to  solve  the 
time-dependent,  3D  Navier-Stokes  equations  for  incompressible  flow. 
Additionally,  the  CFD  model  must  account  for  the  moving  boundary 
corresponding  to  the  parachute  canopy. 
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•  Sophisticated  mesh  generation  algorithms  are  needed  to  automatically  generate 
and  update  the  CFD  mesh  with  moving  boundaries  and  interfaces.  As  the 
parachute  system  deforms  and  travels  through  the  fluid  medium,  the  fluid  mesh 
must  conform  to  the  moving  boundary. 

•  Computational  methods  are  needed  to  simulate  the  effects  of  adverse  wind  fields 
on  parachute  systems.  The  CFD  model  must  also  include  a  general  method  to 
prescribe  wind  conditions  on  the  parachute  model. 

•  Stand-alone  CFD  simulations  are  needed  to  verify  the  fluid  dynamics  model.  To 
verify  the  CFD  model,  stand-alone  simulations  are  performed  with  prescribed 
structural  geometry. 

Fluid-Structure  Interaction 

•  A  robust  method  is  required  to  couple  the  CSD  and  CFD  to  perform  FSI 
simulations.  FSI  simulations  require  the  interchange  of  data  between  the  CSD 
and  CFD  models.  Specifically,  the  CSD  model  must  pass  the  structural 
displacements  and  velocities  to  the  CFD  model  and  the  CFD  model  must  pass  the 
fluid  pressures  to  the  CFD  model.  This  data  interchange  occurs  for  each  time  step 
and  an  iterative  coupling  scheme  is  required  to  insure  convergence  of  the  FSI 
model. 

5.  SUMMARY  OF  RESULTS 

The  goal  of  this  project  has  been  to  develop  computational  tools  for  the  simulation  of 
candidate  airdrop  systems  for  the  New  World  Vistas  (NWV)  Precision  Airdrop  (PAD) 
program.  To  achieve  this  goal,  the  following  results  were  completed: 

•  The  CSD,  CFD,  and  FSI  models  were  formulated  and  implemented  in  a  parallel 
computational  environment. 

•  A  robust  CSD  model  was  formulated  and  implemented.  A  new  wrinkling 
algorithm  was  implemented  in  the  CSD  model.  Details  of  the  wrinkling 
algorithm  and  numerical  results  are  given  in  References  [1]  and  [2]. 

•  A  method  for  modeling  PMAs  was  formulated  and  implemented.  The  PMA  is 
treated  as  a  geometrically  nonlinear  anisotropic  membrane.  Details  of  this 
formulation  and  numerical  results  are  given  in  References  [3]  and  [4]. 

•  Stand-alone  CSD  simulations  representative  of  NWV  PAD  airdrop  systems  were 
performed  to  verify  the  CSD  model.  Numerical  results  from  these  simulations  are 
given  in  References  [5]  and  [6]. 
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•  The  CFD  model  uses  the  Deforming-Spatial-Domain  /  Stabilized  Space  Time 
(DSD/SST)  formulation  developed  by  the  Team  for  Advanced  Flow  Simulation 
and  Modeling  (T*AFSM)  at  Rice  University.  This  approach  is  ideally  suited  for 
parachute  problems  since  they  require  modeling  of  a  moving  boundary  in  the  flow 
field.  Details  of  this  formulation  are  given  in  References  [7]  and  [8]. 

•  A  formulation  for  coupling  the  CSD  and  CFD  models  to  perform  FSI  simulations 
was  developed  and  implemented.  The  formulation  uses  a  tight  coupling  scheme 
that  iterates  on  the  CSD  and  CFD  solutions  at  each  time  step  to  converge  on  the 
coupled  solution.  Details  of  this  formulation  and  simulation  results  are  given  in 
References  [9]  and  [10] . 

•  The  FSI  model  was  validated  by  comparison  of  simulation  results  with  wind 
tunnel  measurements.  Agreement  between  the  simulation  results  and 
measurements  was  very  good.  Results  from  this  validation  study  are  given  in 
Reference  [11]. 

•  FSI  simulations  were  performed  on  candidate  systems  to  evaluate  the  ability  to 
control  these  systems  using  prescribed  riser  length  changes.  The  candidate 
systems  included  both  T-10  and  G-12  canopies.  Results  from  these  simulations 
for  the  T-10  canopy  are  given  in  References  [12]  and  [13]. 

•  A  major  difficulty  was  encountered  in  performing  FSI  simulations  with  control 
operations.  The  control  operations  typically  result  in  large  distortion  of  the 
inflated  canopy  accompanied  by  contact  between  adjacent  gores.  This  difficulty 
is  illustrated  in  Figures  1  and  2. 

•  The  large  distortion  of  the  canopy  and  contact  of  adjacent  gores  generally  resulted 
in  failure  of  the  CFD  model  due  to  gross  distortion  of  the  CFD  mesh.  A  method 
to  overcome  this  difficulty  was  not  developed  during  the  project  duration  and  this 
became  the  most  limiting  factor  for  performing  additional  FSI  simulations. 

•  Our  research  team  will  continue  to  perform  FSI  simulations  of  NWV  PAD 
airdrop  systems  under  a  three-year  DoD  Challenge  Project  titled  “ Airdrop  System 
Modeling  for  the  21st  Century  Airborne  Warrior”. 


5 


VIEW  DEFLECTED  GEOMETRY  (USER  SPECIFIED  FRAME) 


X-AXIS 

Figure  1:  Deformed  Configuration  with  Control  Operation 


Figure  2:  Skirt  Profile  with  Gore  Contact 
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The  dynamic  behavior  of  parachute  systems  is  an  extremely  complex  phenomena  characterized  by  nonlinear, 
time-dependent  coupling  between  the  parachute  and  surrounding  airflow,  large  shape  changes  in  the  parachute,  and 
three-dimensional  unconstrained  motion  of  the  parachute  through  the  fluid  medium.  Because  of  these  complexities, 
the  design  of  parachutes  has  traditionally  been  performed  using  a  semi-empirical  approach.  This  approach  to 
design  is  time  consuming  and  expensive.  The  ability  to  perform  computer  simulations  of  parachute  dynamics 
would  significantly  improve  the  design  process  and  ultimately  reduce  the  cast  of  parachute  system  development 
The  finite  element  formulation  for  a  structural  model  capable  of  simulating  parachute  dynamics  is  presented. 
Explicit  expressions  are  given  for  structural  mass  and  stiffness  matrices  and  internal  and  external  force  vectors. 
Algorithms  for  solution  of  the  nonlinear  dynamic  response  are  also  given.  The  capabilities  of  thestructural  model  are 
demonstrated  by  three  example  problems.  In  these  examples,  the  effect  of  the  surrounding  airflow  is  approximated 
by  prescribing  the  canopy  pressure  and  by  applying  cable  and  payload  drag  forces  on  the  structural  model. 
The  examples  demonstrate  the  ability  to  simulate  three-dimensional  unconstrained  dynamics  beginning  with  an 
unstressed  folded  configuration  corresponding  to  the  parachute  cut  pattern.  The  examples  include  simulations  of 
the  inflation,  terminal  descent,  and  control  phases. 
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Nomenclature 

=  determinant  of  metric  tensors ' A af}  and  Aa/i, 
respectively 

=  covariant  metric  tensors  of  'C  and  °C,  respectively 
=  contravariant  metric  tensors  of  'C  and  °C, 
respectively 

=  material  constitutive  tensor 
=  Cartesian  displacements  of  node  N  on  element 
from  °C  to  'C 

—  Cartesian  velocities  of  node  N  on  clement  in  fC 
=  Cartesian  accelerations  of  node  N  on  element  in  'C 
=  correct  and  estimated  displacement  vector  of  nodes 
from  °C  to  'C,  respectively 
=  Young’s  modulus 

=  unit  base  vectors  in  Cartesian  directions 
=  permutation  symbol 
=  discretized  external  force  on  node  N 
=  polynomial  shape  function  of  xa 
=  thickness  of  membrane 
=  effective  stillness  matrix  for  Newmark  method 
=  estimated  tangent  stiffness  matrix  on  *C 
=  mass  matrix  connecting  nodes  M  and  N 
=  unit  normal  to  'C 

=  Kirchhoff  stress  resultant  tensors  on  'C  and  *C\ 
respectively 

=  net  pressure  on  canopy 
=  Cartesian  components  of  tractions  on  ’C 
=  discretized  internal  restoring  force  on  node  N 
=  time 

=  Cartesian  displacements  of  point  to  'C  from  °C 
=  Cartesian  accelerations  of  point  to  'C  from  °C 
=  nondimensional  curvilinear  coordinates  in  element 
=  Cartesian  coordinates  of  node  N  on  element  in  °C 
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y,  =  Cartesian  coordinates  of  a  point  on  °C 

'  v,*  =  Cartesian  coordinates  of  a  point  on  'C 

c x,  0  =  Newmark  parameters 

'yap  =  Green  strain  tensor 

(#A£>)  =  correction  vector  to  estimate  (*£>} 

A /  =  time  step 

Si,  =  Kronecker  delta 

<5u,  =  virtual  displacements 

r/,  k  =  damping  parameters 

v  =  Poisson’s  ratio 

co  =  damping  ratio  and  natural  frequency  of  a  particular 

mode  of  free  vibration 
p  =  mass  density  of  °C 


I.  Introduction 

PARACHUTE  deployment,  inflation,  and  terminal  descent  are 
extremely  difficult  aerodynamic  phenomena  to  model.  These 
processes  are  governed  by  coupled,  nonlinear,  time-dependent  equa¬ 
tions  for  the  interaction  of  the  fluid  medium  and  parachute.  Be¬ 
cause  of  these  complexities,  it  is  not  surprising  that  the  design  of 
parachutes  has  historically  been  performed  using  semi-empirical 
methods.1 

During  the  last  two  decades,  the  demands  placed  on  parachute  de¬ 
signers  have  increased  tremendously.  Payload  costs  have  increased, 
mission  requirements  have  become  more  stringent,  and  the  flight 
tests  needed  to  design  new  parachute  systems  have  become  more 
costly.  In  light  of  these  demands,  the  traditional  semi-empirical  ap¬ 
proach  to  parachute  design  is  inadequate.  Computational  methods 
have  the  greatest  potential  for  providing  the  necessary  predictive 
models  for  parachute  design. 

The  three  phases  that  occur  during  a  parachute  mission  are  deploy¬ 
ment,  inflation,  and  terminal  descent  with  precision  soft  landing.2 
The  deployment  phase  is  a  highly  nonlinear  process  in  which  the 
parachute  and  suspension  lines  unfold  from  the  deployment  bag.  In 
the  inflation  phase,  airflow  causes  large  canopy  shape  changes  as 
the  parachute  is  stressed  and  begins  to  decelerate.  This  phase  is  also 
highly  nonlinear,  and  maximum  forces  typically  occur  during  this 
time.  The  terminal  descent  phase  corresponds  to  all  behavior  follow¬ 
ing  inflation  and  may  include  various  control  operations  performed 
on  the  parachute  system.  In  general,  all  three  phases  involve  strong 
coupling  between  the  parachute  system  and  surrounding  airflow. 
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Computer  simulations  of  all  three  phases  are  needed  to  fully  char¬ 
acterize  the  behavior  of  a  given  parachute  system  under  given  con¬ 
ditions.  In  general,  comprehensive  simulations  require  coupling  of 
structural  dynamics  (SD)  and  computational  fluid  dynamic  (CFD) 
models  to  capture  the  true  aerodynamic  behavior,  however,  valuable 
information  can  be  obtained  from  either  standalone  SD  or  CFD  sim¬ 
ulations.  The  finite  element  method  (FEM)  is  most  commonly  used 
for  the  SD  model  and  grid-based  CFD  and  vortex  element  meth¬ 
ods  are  both  used  for  fluid  modeling.  Standalone  SD  simulations 
are  performed  by  prescribing  the  canopy  pressure  and  aerodynamic 
drag  on  the  parachute  components  and  are  computationally  efficient 
because  the  SD  model  is  usually  considerably  smaller  than  the  re¬ 
quired  CFD  model.  For  standalone  CFD  simulations,  the  parachute 
geometry  is  prescribed,  and  the  surrounding  airflow  is  calculated. 
The  calculated  surface  loads  are  then  used  to  derive  global  parachute 
motion. 


coordinates  embedded  in  the  reference  surface.  The  nonlinear  dy¬ 
namic  responses  are  integrated  in  time  using  an  implicit  Newmark 
method.19 

In  the  following,  capital  Latin  letters  as  indices  denote  nodal  vari¬ 
ables  on  an  element.  Greek  indices  take  on  the  values  1  and  2  on  the 
surface,  and  lower  case  Latin  indices  take  on  the  values  1,  2,  and 
3.  A  repeated  index  indicates  summation  over  the  range  of  the  in¬ 
dex.  Commas  used  as  subscripts  denote  partial  differentiation  with 
respect  to  the  coordinate  component  with  the  succeeding  index.  A 
special  notation  is  adopted  to  identify  various  configurations  of  the 
membrane:  A  leading  superscript  on  a  variable  denotes  the  surface 
on  which  the  variable  is  defined.  The  reference  surface  where  mag¬ 
nitudes  and  orientations  of  the  variables  are  defined  is  the  initial 
unstrained  state  °C.  Thus,  we  will  arrive  at  a  total  Lagrangian  de¬ 
scription  of  the  nonlinear  equations  at  time  t  on  surface  'C,  or  at 
time  t  +  At  on  *C. 


II.  Previous  Work 

Deployment  has  been  successfully  simulated  by  Purvis3  using  a 
discrete  element  model  similar  to  Sundberg.4  In  that  axisymmet- 
ric  model,  the  partially  deployed  canopy  and  suspension  lines  arc 
modeled  as  elastically  connected,  lumped  mass  nodes,  and  semi- 
empirical  parameters  account  for  aerodynamic  forces.  Although 
improvements  to  this  model  are  possible,  three-dimensional  finite 
element  modeling  of  the  initial  deployment  process  have  not  yet 
been  performed. 

Significant  progress  has  recently  been  made  on  simulation  of 
the  parachute  inflation  and  terminal  descent  phases.  Coupled  infla¬ 
tion  simulations  of  axisymmetric  parachutes  have  been  presented 
by  Haug  et  al.5  and  Stein  et  al.6  Three-dimensional  standalone 
SD  inflation  simulations  have  been  presented  by  Mossecv7  and 
Benney  et  al.8  Standalone  CFD  models  have  been  used  to  predict  the 
flowficld  in  the  terminal  descent  phase  surrounding  a  canopy  with 
prescribed  geometry.  This  type  of  analysis  has  been  performed  for 
round  parachutes,9  round  parachute  clusters,10  and  parafoils.11  • 12 
Standalone  SD  simulations  of  cross  and  round  parachutes  and 
parafoils  in  terminal  descent  have  been  demonstrated  by  Benney 
et  al.8  Coupled  simulations  of  the  terminal  descent  phase  begin¬ 
ning  with  an  inflated  canopy  shape  have  recently  been  performed 
by  Nelsen,13  Stein  et  al.,14  and  Benney  et  al.15 

In  summary,  the  ability  to  perform  standalone  SD  inflation  sim¬ 
ulations  has  been  demonstrated  by  several  researchers.  Standalone 
CFD  predictions  of  the  flowfield  surrounding  a  parachute  with  a  pre¬ 
scribed  rigid  shape  are  generally  available.  Coupled  dynamic  and 
deformable  inflation  simulations  have  been  very  limited.  Coupled 
simulations  of  the  terminal  descent  phase  beginning  with  an  inflated 
canopy  shape  have  recently  been  achieved  by  several  researchers.  To 
date,  detailed  simulation  of  the  initial  deployment  phase  has  not  been 
addressed.  There  still  exists  a  tremendous  need  to  perform  coupled 
dynamic  simulations  of  all  three  phases  of  a  parachute  dynamics 
problem. 

The  purpose  of  this  paper  is  to  present  the  theoretical  foundations 
for  an  SD  model  capable  of  simulating  parachute  deployment,  infla¬ 
tion,  and  terminal  descent.  This  model  is  being  developed  jointly  by 
researchers  at  the  U.S .  Army  Soldier  Systems  Command,  Natick  Re¬ 
search,  Development,  and  Engineering  Center,  and  the  University 
of  Connecticut,  and  a  parallel  version  of  the  model  is  being  cou¬ 
pled  with  a  CFD  model  developed  by  researchers  at  the  Army  High 
Performance  Computing  Research  Center  at  Rice  University.  The 
structural  theory  has  been  incorporated  into  a  finite  element  code 
that  is  highly  tailored  for  parachute  simulations.  Several  standalone 
SD  simulations  are  presented  to  demonstrate  the  capabilities  of  this 
structural  code.  The  example  problems  are  rather  complex,  and  cor¬ 
responding  experimental  data  is  not  readily  available.  The  predic¬ 
tions  of  the  coupled  model  have  recently  been  compared  with  ex¬ 
periments  and  good  agreement  was  obtained.16 

III.  SD  Theory 

In  this  section  we  consider  the  nonlinear  dynamic  behavior  of 
thin,  extremely  flexible,  and  cable-reinforced  membranes  using  the 
FEM.17* 18  Isoparametric  elements19  are  used  to  develop  discretized 
equations  of  motion  admitting  large  dynamic  displacements.  A  total 
Lagrangian19  coordinate  system  is  adopted,  with  curvilinear  local 


A.  Geometric  Interpolation 

The  simultaneous  interpolation  of  geometry  as  well  as  displace¬ 
ments  by  the  same  shape  functions,  that  is,  isoparametric  interpo¬ 
lation,  can  be  used  for  curved  cables  and  membranes.  We  specify 
curvilinear  coordinates  X\  and  on  the  reference  surface  °C  to  be 
nondimensional  natural  coordinates  and  assume  the  Cartesian  co¬ 
ordinates  Vj  and  displacements 'w,  at  a  generic  point  in  the  element 
to  be  given  by 

y,CO  =  HN(xa)YiN  (la) 

'ui{xa)  =  HN(xa)’DiN  (lb) 

where  /f,vCO  are  the  clement  shape  functions.  We  adopt  polyno¬ 
mial  shape  functions  that  map  a  general  quadrilateral  or  triangu¬ 
lar  surface  element  curved  in  three-dimensional  space  onto  two- 
dimensional  squares  or  triangles. 

As  the  surface  of  the  membrane  deforms  from  the  unstrained  state 
°C  to  'C  at  time  t  we  write  the  time-varying  coordinates  as 

'y,-  =  y,-  +  Mi  (2a) 

Because  the  metric  tensor  of  a  differential  arc  length  in  'C  relative 
to  °C  is  given  by20  [using  Eqs.  (I)] 

'  A*fi  =  #>’i.«  *  'yi.fi  —  T*  ’ DiN){YiSf  +  'Dim)  (2b) 

and  because  the  Green  strain  tensor20  of  a  differential  arc  length  on 
'C  relative  to  °C  is 

/  Aap  Aap 
Yafi  -  - 2 - 

jj  u  Yw’DiM  +  Yi^Dis  4-  'Dim'Dim 
=  ttN.ariM.fi - 2 -  W 

the  stress  resultant  tensor  for  a  Hookean  material  is 

'nafi  =  Ca^  (4) 

where  for  an  isotropic  material  in  plane  stress20 
hF"  v  h.F 

-Aal‘Ax'1  +  —  (A^A^  +  A^A^)  (5) 


B.  Equations  of  Motion 

We  will  develop  the  discretized  equations  of  motion  in  terms  of 
an  arbitrary  constitutive  equation  that  will  allow  several  existing 
material  models  to  be  used;  Eq.  (5)  is  a  particular  example.  The 
equations  of  motion  can  be  developed  from  the  principle  of  virtual 
work  as  written  on  the  dynamic  equilibrium  state  'C  but  referred  to 
°C  as17 


I  -  J  J -  'p,Suj)VA  d.t|  dx2 
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where  <5y0/j  is  the  infinitesimal  virtual  strain  corresponding  to  a  small 
virtual  displacement  Su{  superposed  on  'C,  and  the  force  vector 

F\/Md.t|  dx2  =  'p,VA  dxj  d;c2£;  (7) 

has  been  expressed  in  terms  of  Cartesian  components  ' p,*. 

In  a  consistent  mass  approach,  the  accelerations  are  given  by 

'«  -  Hu’Dm  (8) 

Then  substitute  Eqs.  ( 1 ),  (3),  and  (8)  into  Eq.  (6)  to  obtain  for  an  arbi¬ 
trary  virtual  displacement  SDMi  the  following  discretized  equation 
of  motion  for  the  ith  force  component  at  the  Nth  node: 

j  +  'Rm  =  (9a) 

where  the  mass  matrix  is 

Musi}  =  &ijph  If  HNHMVAdxtdx2  (9b) 

the  internal  restoring  force  is 

'Rm  —  j  J  J  riaf*  {Hm.oHn.p  + 

X  (Y\ti  +  Dsti)  (9c) 

and  llie  external  force  is 


C.  Nonlinear  Forced  Vibrations 

Here  we  will  develop  a  solution  scheme  for  predicting  the  non¬ 
linear  dynamic  motions.  The  scheme  is  a  combined  incremental 
and  iterative  method.  We  adopt  a  total  Lagrangian  approach19  in 
which  we  disregard  nonlinear  interactions  of  the  membrane  with 
the  embedding  medium.  If  we  consider  the  surface  *C  at  a  slightly 
later  time  /  -f  A/  from  the  surface  'C  at  time  /  when  we  assume  all 
variables  are  known,  we  can  rewrite  Eq.  (9)  at  time  /  +  Ar  as 

[A/]{#D}  +  (#R}  =  (#F]  (12) 

where  a  bracketed  term,  for  example,  [A/],  denotes  a  matrix,  for  ex¬ 
ample,  Mum}*  and  braces  about  a  term,  for  example,  (*/?),  denotes 
a  vector,  for  example,  *RNi.  The  vectors  (*/?)  and  {*F]  are  obtain¬ 
able  from  Eqs.  (9)  by  replacing  the  superscript '  by #.  Note  that  [ M ] 
is  neither  dependent  on  time  nor  on  displacements  and,  therefore, 
needs  to  be  calculated  only  once. 

If  we  assume  that  the  exact  value  of  [*  D)  at  time  t  4-  Ar  is  given 
by  an  iteration  on  a  prior  estimate  (#D)  to  the  solution,  that  is, 

{#D}  =  (*D)4(#AZ)]  (13a) 

then  we  can  take  Taylor  expansions  in  Eqs.  (12)  about  the  solution 
(#D)  and  obtain  a  Newton-Raphson  iterative  equation22 

{M]\*b)  +  l‘K)(#AD)  =  U'F)  -  {#F}]  (13b) 

in  which  [* F\  —  (*/?)  is  evaluated  at  {* D\  and  (*K)  is  the  tangent 
stiffness  matrix  evaluated  at  (#D) 


' F Si  -  J  I''p,HsVa  d.t,d.i2 


(9d) 


3((^)-{#n) 

K  j  —  - - - 

31  *D) 


(13c) 


Note  that  J  A  in  Eqs.  (9)  can  be  obtained  from  Eq.  (2b)  in  terms  of 
derivatives  of  the  shape  functions  and  the  nodal  coordinates,  thus, 
further  complicating  the  integration  required  to  form  the  element 
matrices. 

Equations  (9)  represent  a  set  of  ordinary  differential  equations 
of  dynamic  equilibrium  at  time  t  of  every  degree  of  freedom  on 
'C.  Equations  (9)  are  nonlinear  in  that  1)  material  nonlinearities  are 
possible  in  the  restoring  force  because  'nafi  can  be  replaced  by 
any  constitutive  equation;  2)  1  RNi  is  nonlinear  even  if  the  consti¬ 
tutive  equation  is  linear  because  includes  finite  displacement 
components  '£>>,,•;  and  3)  nonconservative  loads  arc  possible  in  the 
external  force ’  FS;  if  the  magnitude  or  direction  of  'p,-  is  dependent 
on  displacements,  for  example,  pressure. 

In  pressurized  membranes,  nonconservative  loads  occur  during 
large  displacements.  The  tractions  'p,-  due  to  net  pressure  p  normal 
to  the  membrane  surface  'C  are 

'p,  =  ■J'A/Ap(Nei)  (10) 

where  the  dot  between  'N  and  e,-  represents  the  inner  vector  product 
of  the  unitjiormal  to  the  surface  with  the  Cartesian  base  vector  e,-. 
Because  'N  is  defined  in  terms  of  the  cross  product  of  the  deformed 
base  vectors,  we  replace  'yf  =  y}  by  Eqs.  (1)  and  find  that  for 
pressure  loading  Eq.  (9d)  becomes 


If  conservative  loads  are  considered,  d[* F}!d{* D)  =  0.  If  non¬ 
conservative  loads,  such  as  pressure,  are  considered,  {*  F)  is  obtained 
from  Eq.  ( 1 1 )  and,  therefore, 


a  l*h 


=  Peiji  (YQk  +  *DQk)j J  FIs  Hq.\ Hm.z  d.ti  d.r2 


af  *D) 

+  (YQk-,DQk)  J j  HnH„aHq.2  djt,djf2 


(14a) 


Because  of  the  permutation  symbol  ey,-,  we  see  that  Eq.  (14a)  will 
contribute  nonsymmetric  terms  to  the  tangent  stiffness  matrix.  It 
is  often  convenient  to  achieve  a  symmetric  stiffness  matrix  and, 
therefore,  in  one  solver  described  in  Sec.  IV,  the  nonsymmetric 
terms  in  Eq.  (14a)  were  disregarded,  and  instead  the  load  {#F}  in 
Eq.  ( 1 2)  was  treated  iteratively  to  partially  account  for  that  nonlinear 
effect.  The  full  nonsymmetric  terms  are  included  in  another  solver 
described  in  Sec.  IV. 

Thus,  the  principal  term  in  the  tangent  stiffness  matrix  is 


t'K] 


3f*) 
3  [*£>) 


(14b) 


and  is  obtained  by  taking  the  indicated  derivatives  of  Eq.  (14b), 


’ Fw  =  ptkjiyj  j  HnHqj  HMtldx\  dxij 

x  {YQk  Y\h  4-  YQk'DMj  4*  YQk'DMi  +' DQk' D^j)  (11) 

where  ek)i  is  the  permutation  symbol  introduced  because  of  the  cross 
product.  We  see  that  'p/  and,  hence, ' Fs;  are  indeed  nonlinear  in  the 
displacement  components  'Dm;. 

To  evaluate  the  matrices  in  Eqs.  (9),  we  must  integrate  products 
of  the  shape  functions  HM{. O  and  their  derivatives  over  the  area  of 
each  element.  Because  of  the  complexities  of  fhe  those  products,  the 
integrations  must  be  done  numerically.21  Gaussian  quadrature  can 
be  used  for  the  quadrilateral  elements,  and  Gauss-Hammer  formulas 
can  be  used  for  the  triangular  elements. 


KMmi  =  fj  +  HN.aHM.f) 

+  X-C^^HR.aHN.„  +  H„.aHR,l>)(Yl'l+'DRi)(HQ,nHM,( 

+  +  *DQj)  Jyjd.T,  d.r2  (15) 

D.  Newmark's  Method19 

Given  the  solution  on  'C  at  time  t  and  the  linearized  iterative  form 
of  the  equations  of  motion,  Eqs.  (13-15),  we  seek  an  incremental 
solution  for  the  iteration  {*  AZ))  to  the  solution  on  "C,  Toward  that 
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end,  we  assume  distributions  for  the  velocities  {#D}  and  accelera¬ 
tions  (*D)  during  the  time  interval  from  t  to  t  +  A/  such  that 

(*Z>)  =  {'/))  +  [(1  -«){'/)) +*{'/))]  A/  (16a) 

[#D]  =  {'D)  +  f£>}Ar  +  [{{  -  {*){'£>)  +  fil'd)] A/2  (16b) 

where  a  and  p  are  parameters  dictating  the  assumed  distribution  of 
acceleration  over  the  interval;  different  explicit  and  implicit  formu¬ 
lations  can  be  developed  with  different  parameters.  If  a  is  equal  to 
or  greater  than  0.5  and  if  p  is  equal  to  or  greater  than  0.25 (0.5  + a)2, 
the  formulation  is  implicit.19 

In  the  implicit  formulation  Eq.  ( 1 6b)  is  solved  for  {*£)}  in  terms  of 
values  at  time  t  and  in  terms  of  {*D).  On  substitution  into  Eqs.  (13), 
an  equivalent  static  problem  is  derived  at  each  time  step  in  the  form 

[KW]l#AD)  =  {£)  —  {/?efrl  07a) 

where 

[ATe(I]  =  ('An  +  [W]/(M»2)  07b) 

is  the  effective  stiffness  matrix  and  is  the  similarly 

modified  load  vector. 

E.  Damping 

Artificial  damping  of  results,  particularly  of  higher  modes  of  vi¬ 
bration,  can  be  realized  by  varying  the  Newmark  parameter  or  in 
Eqs.  (16).  If  a  =0.5.  no  artificial  damping  is  introduced.  For  in¬ 
creasingly  larger  values  of  or  beyond  0.5,  increasingly  rapid  artificial 
damping  is  achieved.  This  would  be  desired  if  steady-state,  rather 
than  transient,  solutions  are  sought. 

If  real  damping  is  present  in  the  system,  or  an  alternative  artificial 
damping  mechanism  is  sought,  Rayleigh  damping  proportional  to 
the  stiffness  and  mass  matrices  can  be  introduced.  Add  a  damping 
term  [*C](#Z))  to  Eq.  (15)  with  the  damping  matrix  assumedto  have 
the  form23 

[#C]  =  +  (18a) 

and  then  add  appropriate  terms  to  [/fcfr]  and  (/?crr|.  For  larger  values 
of  the  lower  transient  modes  are  more  heavily  damped.  For  larger 
values  of  X,  the  higher  transient  modes  are  more  heavily  damped. 
The  damping  ratio  £  relative  to  critical  damping  of  a  particular  modal 
frequency  co  is  given  in  terms  of  tj  and  X  by19 

£  =  ri/(  2o))  +  Xw/2  (18b) 

IV.  Computer  Implementation 

The  SD  theory  described  in  the  preceding  section  has  been  im¬ 
plemented  in  a  general  purpose  FORTRAN  computer  code,  called 
TENSION,  for  analysis  of  the  nonlinear  dynamic  response  of  ca¬ 
ble  and  membrane  structures.  Numerous  special  features  have  been 
added  to  this  code  to  facilitate  computer  simulation  of  parachute  dy¬ 
namics.  A  custom  pre-  and  postprocessor  written  in  MATLAB®  has 
also  been  developed  to  provide  for  state-of-the-art  model  generation 
and  visualization  of  the  simulation  results. 

TENSION  currently  has  three  dynamic  solvers.  Two  implicit 
solvers  are  based  on  Newmark’s  method,  and  the  third  is  based 
on  the  explicit  central  difference  method,19  which  is  a  conditionally 
stable  method.  The  difference  between  the  two  implicit  solvers  is  in 
the  method  used  to  solve  the  algebraic  equations  given  by  Eq.  (17a). 
The  first  uses  a  direct  factorization  method21  with  a  profile  reduc¬ 
tion  algorithm24  and  is  only  applicable  to  symmetric  systems,  and 
the  second  employs  an  iterative  method25  26  that  can  be  used  for 
nonsymmetric  systems  of  equations. 

V.  Results 

In  this  section,  three  example  problems  are  presented  to  de¬ 
monstrate  the  capabilities  of  TENSION  for  simulating  parachute 
dynamics.  All  three  problems  are  performed  as  standalone  SD  sim¬ 
ulations.  The  fluid  effects  are  approximated  by  prescribing  rep¬ 
resentative  pressures  on  the  canopy  membranes  and  by  including 
velocity  dependent  fluid  drag  forces  on  cables  and  concentrated 
masses.8  All  examples  utilize  a  linear  elastic  material  model  with 


b) 

Fig.  1  Initial  folded  configuration  of  cross  canopy,  a)  folding  of  one 
arm  (upside  down)  and  b)  entire  canopy. 

constant  material  properties.  The  properties  used  in  these  sim¬ 
ulations  are  £  =  4.32x  106  lb/ft2  =0.207  GPa,  v  =  0.3,  mem¬ 
brane  thickness  =  0.03  mm  (0.0001  ft),  cable  cross-sectional  area  = 
9.29  mm2  (0.0001  ft2),  fabric  density  =  309.3  kg/m3  (0.6  slugs/ft3), 
air  density  =  1.222  kg/m3  (0.00237  slugs/ft3),  and  gravitational  ac¬ 
celeration  =  9.81  m/s2  (32.2  ft/s2). 

A.  Opening  and  Control  of  a  Cross-Type  Canopy 

This  example  demonstrates  the  ability  of  TENSION  to  simulate 
the  opening  of  a  full  three-dimensional  unconstrained  model  of  a 
cross-type  parachute  in  an  initially  folded  configuration  and  subse¬ 
quent  control  of  the  parachute  motion  through  user-defined  changes 
of  the  control  line  lengths.  The  four  arms  of  the  cross  canopy  are 
initially  folded  over  the  canopy  center,  as  shown  for  one  arm  in 
Fig.  1.  In  this  simulation,  the  model  is  subjected  to  constant  canopy 
pressure,  cable  and  payload  drag  that  is  dependent  on  the  square  of 
the  relative  velocity,  and  gravity.  Time-dependent  mass  proportional 
damping  is  applied  during  the  initial  opening  to  stabilize  the  sim¬ 
ulation,  then  reduced  after  partial  opening  is  achieved.  The  initial 
geometry  and  deformed  geometries  at  three  times  during  opening 
are  shown  in  Fig.  2. 

Following  the  opening,  the  model  approaches  a  steady-state  de¬ 
scent.  Next,  four  separate  control  line  length  changes  are  prescribed 
to  simulate  the  parachute  flying  in  a  box  pattern.  Control  line  changes 
are  effected  in  the  FEM  model  by  prescribing  time-dependent  un¬ 
stressed  cable  lengths.  Each  control  operation  consists  of  lengthen¬ 
ing  a  control  line  for  10  s,  then  retracting  it  for  10  s  prior  to  initiating 
the  next  control  operation.  The  .r,  y,  and  z  velocities  of  the  payload 
for  the  entire  simulation  are  shown  in  Fig.  3.  Prior  to  control,  the 
horizontal  (.t,  y)  velocities  are  zero,  and  terminal  velocity  is  ap¬ 
proached  in  the  vertical  z  direction.  Control  operations  begin  at  20  s 
and  induce  the  horizontal  velocities  shown.  The  trajectory  of  the 
payload  is  shown  in  Fig.  4.  Here,  the  initial  descent  prior  to  control 
and  the  four  periods  of  control  are  clearly  evident.  The  deformed 
shape  of  the  parachute  model  initiating  the  first  control  operation  is 
shown  in  Fig.  5. 

B.  Opening  of  a  Ram-Air  Parafoil  System 

This  model  simulates  a  ram-air  parafoil  (wing)  that  is  similar 
to  that  flown  under  the  U.S.  Army  GPADS27  and  NASA  X-3828 
programs.  For  simplicity  and  reduction  of  computational  time,  the 
parafoil  is  modeled  in  half-plane  symmetry.  The  wing  is  initially  in 
a  flat  configuration  that  represents  the  cut  pattern  geometry.  Figure  6 
shows  the  full  three-dimensional  initial  geometry. 

The  pressure  distribution  used  in  this  simulation  is  constant  for 
all  time.  This  is  possible  due  to  time-dependent  mass  proportional 
damping  that  is  slowly  minimized  during  the  early  portion  of  the 
simulation.  The  top,  bottom,  and  side  faces  are  supplied  different 
pressure  distributions.  The  system  is  acted  on  by  a  gravitational 
load  in  the  negative  z  direction,  mass  proportional  damping,  and 
velocity-dependent  cable  drag. 

The  time-dependent  mass  proportional  damping  is  ramped  from 
an  initially  high  value  to  a  low  value  over  the  first  5  s  of  the  simu¬ 
lation.  These  values  allow  for  a  slow  oozing  of  the  parafoil  into  the 
inflated  shape  and  a  method  of  allowing  the  suspension  lines,  all  of 
which  start  in  a  slack  state,  to  become  axially  loaded.  The  program 
allows  the  user  to  define  the  cable  length  between  two  nodes  to  be 
longer  than  the  actual  distance  between  the  nodes  and  to  prescribe 
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Fig.  2  Initial  configuration  and  deformed  shapes  during  initial  opening  of  cross  parachute. 
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Fig.  3  Payload  velocity  of  cross  parachute. 


Fig.  4  Trajectory  of  payload  of  cross  parachute. 


that  the  cable  is  in  a  slack  state.  This  is  a  critical  capability  that  al¬ 
lows  the  user  to  run  the  same  model  (same  nodal  coordinates)  with 
different  sets  of  suspension  line  and  bridle  lengths.  The  user  can 
also  prescribe  changes  in  cable  lengths  in  time  that  could  be  used  to 
predict  the  effect  on  the  systems  performance  to  changes  in  rigging 
angle.  As  already  mentioned,  all  suspension  lines  are  initially  in  a 
slack  state  with  a  user-defined  cut  length.  This  also  assists  in  model 
development  due  to  the  relatively  simple  geometry  and  the  ability  to 
incorporate  actual  cut  patterns  for  the  canopy.  The  time -dependent 
mass  proportional  damping  allows  for  the  smooth  transition  to  an 
inflated  shape  and  transition  into  steady-state  glide.  Figure  7  shows 
a  sequence  of  snapshots  during  the  initial  inflation  process. 


The  system  slowly  reaches  a  steady-state  glide  after  approxi¬ 
mately  40  s.  Figure  8  shows  the  y  and  z  positions  of  the  payload 
and  node  A  vs  time.  At  a  simulation  time  of  45  s,  six  trailing-edge 
control  lines  are  shortened  to  a  predefined  length  over  a  5-s  time 
interval  to  simulate  a  flare  manuver.  The  pressure  loading  on  the 
system  during  a  flare  will  most  definitely  change  over  time;  how¬ 
ever,  for  this  simulation,  the  pressure  is  assumed  constant.  The  y 
and  z  velocities  of  the  payload  and  node  A  are  shown  in  Fig.  9, 
which  clearly  shows  the  effect  of  the  flare  maneuver. 

This  example  demonstrates  the  ability  of  TENSION  to  simulate 
the  inflation  of  a  parafoil  starting  from  cut  patterns  to  a  steady-state 
glide  with  the  potential  of  the  user  modifying  the  rigging  angle  in 
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Fig.  5  Deformed  shape  of  cross  parachute  during  first  control  operation. 


Fig.  6  Initial  geometry  of  parafoil  model.  Fig.  8  Position  vs  time  curves  for  parafoil  simulation. 
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Fig.  7  Deformed  shapes  during  parafoil  inflation. 


time  to  optimize  design  and  a  flare-type  control  maneuver.  Gener¬ 
alization  to  a  fully  three-dimensional  unconstrained  model  for  turn 
performance  predictions  is  expected  to  be  a  small  modification  that 
would  simply  increase  the  overall  model  size. 

C.  Opening  of  a  Round  Canopy 

The  third  example  is  a  generic  round  ring-slot  canopy.  The  model 
consists  of  28  gores  and  6  ring  slots.  The  model  starts  in  a  folded 
configuration,  as  shown  in  Fig.  10.  The  simulation  begins  with 
an  initial  horizontal  velocity  on  ail  node  points  in  the  negative 
x  direction.  Note  that  all  28  gores  are  folded  along  the  gore  cen¬ 
terlines  in  a  configuration  that  is  close  to  that  obtained  at  line 
stretch.  The  model  is  fully  three-dimensional  and  unconstrained. 
The  model  is  subjected  to  a  vertical  gravitational  loading,  user- 
defined  time-dependent  pressure  distribution,  velocity-dependent 
cable  drag,  and  velocity-dependent  payload  drag.  Time-dependent 
mass  proportional  damping  is  minimally  present  over  the  first  4  s  of 
the  simulation. 

A  time-dependent  outward  acting  pressure  distribution  is  si¬ 
multaneously  applied  to  all  membrane  elements.  Gravity  is  acting 
throughout  the  simulation.  The  seam  cables,  suspension  line  cables. 
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Fig.  9  Velocity  vs  time  curves  for  parafoil  simulation. 


Fig.  10  Initial  folded  configuration  of  ring-slot  parachute. 
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Fig.  1 1  Deformed  shapes  of  ring-slot  parachute  during  inflation  (0.03- 
0.48  s). 


Fig.  12  Deformed  shapes  of  ring-slot  parachute  during  inflation  (0.51- 
0.96  s). 


Fig.  13  Velocity  vs  time  curves  for  ring-slot  parachute. 


Fig.  14  Final  shape  of  ring-slot  parachute. 


and  payload  are  subjected  to  velocity-dependent  drag  loads  for  all 
times,  which  assists  in  the  trajectory  changing  from  pure  horizontal 
motion  to  pure  vertical  motion. 

A  sequence  of  payload- fixed  snapshots  are  shown  in 
Figs.  11  and  12  for  time  intervals  of  0.03-0.48  and  0.51-0.96  s, 
respectively.  These  snapshots  clearly  show  the  three-dimensional 
capabilities  of  TENSION.  For  example,  the  uppermost  suspension 
lines  become  stressed  and  elongated  during  the  inflation  before  the 


lowermost  suspension  lines.  The  lower  suspension  lines  also  appear 
to  be  most  visibly  affected  by  the  fluid  drag.  Figure  13  shows  the 
horizontal  and  vertical  velocities  of  the  payload  and  node  A  vs  time. 
The  system’s  final  shape  is  shown  in  Fig.  14. 

This  model  demonstrates  the  ability  of  TENSION  to  simulate  a 
round  canopy  with  an  initial  horizontal  velocity  and  folded  config¬ 
uration  through  opening  and  ultimately  into  a  steady-state  vertical 
terminal  configuration. 
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VI.  Conclusions 

The  finite  element  formulation  for  an  SD  model  capable  of  simu¬ 
lating  the  nonlinear  dynamic  behavior  of  highly  flexible  structures 
comprising  membranes  and  cables  has  been  presented.  This  model 
is  currently  being  used  to  simulate  the  inflation,  terminal  descent, 
and  control  of  a  variety  of  parachute  systems.  Three  example  prob¬ 
lems  were  presented  to  demonstrate  the  capabilities  of  this  model 
for  simulating  parachute  dynamics.  In  these  examples,  the  effect  of 
the  surrounding  airflow  was  approximated  by  prescribing  the  mem¬ 
brane  pressure  and  cable  and  payload  drag  on  the  structural  model. 
Coupling  of  the  structural  model  with  CFD  programs  has  been  pre¬ 
sented  in  other  publications.14,15  The  ultimate  goal  of  these  efforts 
is  to  be  able  to  design  and  optimize  parachute  systems  using  com¬ 
puter  simulations,  reduce  the  life  cycle  costs  of  airdrop  systems,  and 
provide  an  airdrop  virtual  proving  ground. 
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ABSTRACT 

We  present  a  strategy  for  carrying  out  3-D  simu¬ 
lations  of  parachute  fluid-structure  interaction,  and 
demonstrate  the  strategy  for  simulations  of  airdrop 
performance  and  control  phenomena  in  terminal  de¬ 
scent.  The  strategy  uses  a  stabilized  space-time  for¬ 
mulation  of  the  time-dependent,  3-D  Navier-Stokes 
equations  of  incompressible  flows  for  the  fluid  dy¬ 
namics  solution.  A  finite  element  formulation  de¬ 
rived  from  the  principle  of  virtual  work  is  used  for 
the  parachute  structural  dynamics.  Coupling  of  the 
fluid  dynamics  with  the  structural  dynamics  is  im¬ 
plemented  over  the  the  fluid-structure  interface, 
which  is  the  parachute  canopy  surface.  Large  de¬ 
formations  of  the  structure  are  handled  in  the  fluid 
dynamics  mesh  using  an  automatic  mesh  moving 
scheme.  The  strategy  is  demonstrated  to  simulate 
the  terminal  descent  behavior  of  a  standard  US  Army 
personnel  parachute  system.  Also,  preliminary  sim¬ 
ulations  are  presented  for  the  response  of  the 
parachute  to  a  variety  of  “riser  slips,”  which  can 
provide  the  parachute  system  with  limited  maneu¬ 
verability. 

INTRODUCTION 

Airdrop  technology  is  a  vital  Department  of  De¬ 
fense  (DoD)  capability  to  the  rapid  deployment  of 
warfighters,  ammunition,  equipment  and  supplies. 

^This  paper  is  declared  a  work  of  the  U.S.  Government  and 
is  not  subject  to  copyright  protection  in  the  United  States. 


In  addition,  airdrop  of  food,  medical  supplies,  and 
shelters  for  humanitarian  relief  efforts  are  increasing 
in  demand.  A  team  of  researchers  centered  at  the 
US  Army  Soldier  and  Biological  Chemical  Command 
(SBCCOM),  Natick  Soldier  Center  (Natick)  are  ac¬ 
tively  involved  in  developing  new  technologies  to  ad¬ 
vance  DoD  airdrop  capabilities.  An  important  tech¬ 
nology  being  developed  is  the  numerical  modeling 
of  parachutes  and  airdrop  system  performance  uti¬ 
lizing  high  performance  computing  (HPC)  resources. 
Parachutes  and  airdrop  systems  have  been  tradition¬ 
ally  developed  by  time  consuming  and  costly  full- 
scale  testing.  The  capability  of  using  HPC  tools 
to  model  and  develop  parachutes  and  airdrop  sys¬ 
tems  will  greatly  reduce  the  time  and  cost  of  full- 
scale  testing,  help  to  minimize  the  life  cycle  costs  of 
airdrop  systems,  assist  in  the  optimization  of  new 
airdrop  capabilities,  and  provide  an  airdrop  virtual 
proving  ground  environment. 

These  modeling  efforts  generally  involve  the  cou¬ 
pling  of  computational  fluid  dynamics  (CFD)  codes 
and  structural  dynamics  (SD)  codes.  These  coupled 
fluid-structure  interaction  (FSI)  models  are  required 
to  capture  the  physics  of  highly  complex  parachute 
dynamic  phenomena.  Applications  and  focus  areas 
include  most  airdrop  systems  and  associated  phe¬ 
nomena.  This  paper  will  focus  on  the  US  Army’s 
mass  assault  personnel  parachute  in  terminal  de¬ 
scent  and  with  riser  control  inputs.  The  T-10 
parachute  system  consists  of  a  flat  extended  skirt 
canopy,  30  suspension  lines  and  four  3-foot  risers 
attaching  to  the  paratroopers  harness.  A  detailed 
description  of  the  structural  model  of  the  canopy 
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and  the  associated  fluid  model  for  the  surrounding 
airfield  will  be  presented. 

Paratroopers  are  trained  at  the  basic  airborne 
school  at  Fort  Benning,  Georgia.  The  three  week 
course  includes  training  in  controlling  a  T-10 
parachute  (the  current  mass  assault  parachute  for 
US  Army  troops)  by  use  of  riser  slips.  Riser  slips 
involve  physically  reaching  up  and  grabbing  a  riser 
and  pulling  it  down  to  modify  the  canopy  shape  and 
induce  glide.  Riser  slips  are  used  to  avoid  close  con¬ 
tact  with  other  jumpers,  avoid  ground  obstacles,  and 
to  “pull  a  two  riser  slip  into  the  wind”  and  hold  it 
approximately  150  to  200  feet  above  the  ground  to 
minimize  relative  horizontal  ground  impact  speeds. 
Experienced  paratroops  can  control  the  T-10  that 
is  typically  a  non-controllable  round  parachute. 

Riser  slips  are  used  routinely  every  day  but  are 
not  fully  understood.  The  effects  of  riser  slip  length, 
jumper  weight,  and  other  variables  on  the  relative 
horizontal  velocity  of  the  system  are  not  well  known 
beyond  jumpers  practical  experience.  The  US  Army 
and  US  Air  Force  are  exploring  the  use  of  riser  slips 
for  use  in  high-altitude  resupply  cargo  systems  to 
allow  for  precision  landing  (i.e.,  target  accuracy)  for 
various  weight  payloads.  These  are  autonomously 
controlled  round  parachutes  in  which  the  risers  are 
replaced  by  actuators  and  which  utilize  GPS  as  the 
primary  navigation  sensor.1  These  systems  require 
knowledge  of  achievable  performance  in  order  for 
control  optimization  strategies  to  be  developed.  Cur¬ 
rently,  the  only  way  to  determine  the  performance 
characteristics  is  through  extensive  instrumented 
testing  of  full-scale  systems. 

In  addition,  the  Army  is  exploring  an  Advanced 
Tactical  Parachute  System  (ATPS)  to  replace  the 
T-10.  The  canopy  chosen  will  not  use  a  T-10  and 
therefore,  the  systems  response  to  a  riser  slip  will 
need  to  be  determined  and  measured  for  compar¬ 
isons  to  the  current  T-10  to  assess  potential  opera¬ 
tional  impacts.  Once  again,  this  will  require  exten¬ 
sive  testing.  The  FSI  capabilities  being  developed 
are  being  applied  to  these  two  real  and  current  ap¬ 
plications.  Concurrently,  and  reported  in  this  paper, 
the  prediction  of  the  performance  of  a  T-10  in  termi¬ 
nal  descent  and  during  a  riser  slip  will  be  presented 
as  a  base  line  example.  A  series  of  detailed  mea¬ 
surements  of  T-10  parachutes  with  personnel  pulling 
riser  slips  is  being  conducted  in  partnership  with  the 
US  Army  Yuma  Proving  Ground  (the  primary  test¬ 
ing  agency  for  such  systems).  However,  the  data  was 
not  available  when  this  paper  was  written. 

The  ultimate  goal  of  this  work  is  to  utilize  these 
tools  to  minimize  feasibility  testing  of  new  systems, 


assist  in  determining  the  capabilities  of  all  systems, 
and  ultimately  to  reduce  the  life  cycle  cost  of  all 
DoD  Airdrop  Systems.  This  paper  will  provide  an 
overview  of  the  FSI  theory  and  finite  element  strat¬ 
egy  being  developed  and  provide  a  status  report  on 
the  ability  to  predict  T-10  performance  in  terminal 
descent  and  during  controlled  riser  deformations. 

MODELING  STRATEGY 
Fluid  Dynamics 

For  the  fluid  dynamics,  the  flow  is  assumed  to 
be  at  a  low  speed  and  thus  the  Navier-Stokes  equa¬ 
tions  of  incompressible  flows  are  utilized.  To  han¬ 
dle  the  time-variant  spatial  domains  encountered  in 
parachute  problems,  we  employ  the  Deforming- 
Spatial-Domain/Stabilized  Space-Time  (DSD/SST) 
formulation2, 3  finite  element  formulation.  In  this 
formulation,  the  finite  element  interpolation  polyno¬ 
mials  are  functions  of  both  space  and  time,  and  the 
stabilized  variational  formulation  of  the  problem  is 
written  over  the  associated  space-time  domain.  This 
stabilized  formulation  automatically  takes  into  ac¬ 
count  deformations  in  the  spatial  domain  and  pro¬ 
tects  the  computation  against  numerical  oscillations. 
This  method  has  been  applied  to  a  large  number  of 
problems  with  moving  boundaries  and  interfaces. 

Structural  Dynamics 

The  parachute  structure  is  represented  as  a  “ten¬ 
sion  structure”  composed  of  membranes,  cables,  and 
point  masses.  For  our  problems,  the  structure  un¬ 
dergoes  large  displacements  and  large  rotations,  but 
small  strains.  Thus,  constitutive  relationships  are 
used  assuming  the  materials  are  Hookean  with  linear- 
elastic  properties.  A  finite  element  formulation  de¬ 
rived  from  the  principle  of  virtual  work  is  used  for 
the  structural  dynamics  (SD).4,5  Finite  displace¬ 
ments  of  the  structure  are  taken  into  account  by 
using  a  total  Lagrangian  description  of  the  prob¬ 
lem.  In  addition  to  membrane  and  cable  elements, 
a  variety  of  parachute-specific  features  have  been  in¬ 
corporated  into  the  SD  solver  to  include  wrinkling 
models6  (which  eliminate  compressive  stresses)  and 
time-variant  cable  lengths7  (for  modeling  of  control 
line  pulls,  “riser-slips”,  reefing,  etc.). 

Mesh-Moving  Strategy 

We  use  an  automatic  mesh-moving  scheme  to 
handle  changes  in  the  spatial  domain  due  to 
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parachute  canopy  deformations.  In  this  scheme,  the 
fluid  mesh  is  treated  as  a  linearly-elastic  “pseudo¬ 
solid”  that  deforms  as  dictated  by  the  motion  of 
the  surface  boundaries  of  the  fluid  domain.8  This 
scheme  introduces  an  additional  computational  cost 
associated  with  the  mesh-moving  equations,  but  is 
well  suited  for  handling  the  complex  geometries  and 
arbitrary  motions  for  this  class  of  problems. 

Fluid-Structure  Coupling 

The  fluid-structure  coupling  occurs  at  the  FSI 
interface,  which  is  in  this  case  the  parachute  canopy 
surface.  We  use  an  iterative  coupling  approach,  with 
individual  systems  of  equations  being  solved  for  the 
fluid  and  the  structure.  Coupling  is  achieved  through 
the  transfer  of  FSI  information  between  the  fluid  and 
structure  within  an  iteration  loop,  with  multiple  it¬ 
erations  improving  the  convergence  of  the  coupled 
system.  Displacements  and  displacement  rates  from 
the  SD  solution  are  treated  as  boundary  conditions 
in  the  mesh  moving  scheme  and  CFD  solver,  respec¬ 
tively.  In  return,  parachute  surface  tractions  from 
the  fluid  are  used  as  distributed  forces  in  the  SD 
solver.  For  the  applications  presented  in  this  paper, 
we  transfer  only  the  pressure  contribution  from  the 
CFD  solution  to  the  SD  solver. 

In  the  example  problems,  we  use  different  meshes 
to  represent  the  parachute  canopy  in  the  SD  and 
CFD  models.  Biquadradic  elements  define  the 
canopy  in  the  SD  mesh  (to  effectively  represent  cur¬ 
vatures  in  the  membranes),  whereas  triangular  ele¬ 
ments  define  the  canopy  surface  in  the  CFD  mesh 
(in  order  to  utilize  automatic  mesh  generation  soft¬ 
ware).  Coupling  information  is  transferred  between 
the  incompatible  meshes  using  a  least-squares  pro¬ 
jection  strategy.  Figure  1  depicts  the  incompatible 
meshes  used  for  the  example  problems  to  be  pre¬ 
sented.  Here,  the  parachute  canopy  is  represented 
with  9-noded  membranes  elements  in  the  SD  mesh 
and  with  3-noded  triangles  in  the  CFD  mesh. 

T— 10  PARACHUTE  APPLICATIONS 

The  described  FSI  modeling  strategy  is  being 
demonstrated  for  a  variety  of  parachute  applications 
to  include  the  the  terminal  descent  of  the  Army’s  T- 
10  personnel  parachute  system.9-13  The  T-10  is  a 
“flat  extended  skirt  canopy”  composed  of  a 
35-foot  diameter  ( Dc )  canopy  and  30  suspension  lines 
each  29.4  feet  long.  The  canopy  is  called  a  “flat  ex¬ 
tended  skirt  canopy”  because  in  its  constructed  (or 
nonstressed)  configuration  it  is  composed  of  a  main 


CFD  mesh  SD  mesh 


Figure  1.  Incompatible  meshes  for  T-10  parachute. 

circular  section  with  a  circular  vent  at  the  apex  and 
an  inverted  flat  ring  section,  which  lies  under  the 
main  section  and  is  connected  to  the  main  section 
at  the  outer  radius.  The  lines  are  connected  to  the 
payload  (i.e.,  paratrooper)  with  four  risers.  The  sus¬ 
pension  lines  continue  as  30  gore-to-gore  reinforce¬ 
ments  through  the  parachute  canopy  and  meet  at 
the  apex.  For  the  T-10,  the  vent  diameter  is  0.1  Dc 
and  the  width  of  the  skirt  is  also  0ADc.  Figure  2 
shows  a  “blown-out”  view  for  the  SD  mesh. 


Figure  2.  SD  mesh  for  T-10  parachute. 
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T-10  Parachute  in  Terminal  Descent 

The  following  example  is  for  the  FSI  of  a  Army 
personnel  parachute  in  terminal  descent  freefall.  The 
simulation  process  involves  three  main  steps. 

Stand-alone  structural  dynamics  simulation : 

The  SD  model  is  broken  into  six  distinct  material 
groups;  one  membrane  group,  three  cable  groups,  a 
truss  group,  and  a  concentrated  mass  group.  The 
membrane  group  defines  the  parachute  canopy.  We 
have  distinct  cable  groups  for  the  suspension  lines, 
the  canopy  radial  reinforcements,  and  the  risers.  The 
truss  and  concentrated  mass  groups  define  the  pay- 
load.  The  SD  mesh  consists  of  3,593  nodes,  780  nine- 
noded  (i.e.,  biquadratic)  membrane  elements  for  the 
canopy  surface,  and  1,196  tw'o-noded  cable  elements, 
and  4  concentrated  mass  points.  The  parachute  sys¬ 
tem  is  represented  by  linearly  elastic  materials,  with 
properties  representative  of  a  T-10.  A  stand-alone 
damped  dynamic  simulation  is  conducted  for  the 
T-10  parachute  model  to  inflate  the  canopy  under 
a  prescribed  differential  pressure  of  0.5  lb/ft2.  For 
the  stand-alone  simulation,  the  four  payload  node 
points  are  fixed  in  space  and  all  other  nodes  are  un¬ 
constrained.  The  equilibrium  shape  for  the  inflated 
parachute  is  depicted  in  Figure  1  (right).  This  equi¬ 
librium  solution  is  used  as  the  initial  condition  for 
the  SD  solver  in  the  subsequent  FSI  simulation. 

Stand-alone  fluid  dynamics  simulation : 

A  3-D  mesh  with  tetrahedral  elements  is  gen¬ 
erated  (with  149,253  nodes  and  888,344  elements) 
using  the  inflated  canopy  from  the  stand-alone  SD 
simulation  as  an  interior  boundary.  Initial  unsteady 
flow  solutions  are  obtained  for  the  fixed  canopy  con¬ 
figuration  at  a  Reynolds  number  of  5  x  106  using  a 
stabilized  semi-discrete  formulation.14  For  the  flow 
simulations,  the  parachute  canopy  is  treated  as  a 
zero-porosity  material  and  the  parachute  canopy  is 
assigned  a  no-slip  boundary  condition.  The  inflow 
boundary  below  the  parachute  is  assigned  a  pre¬ 
scribed  velocity  condition  equivalent  to  the  desired 
Reynolds  number.  Side  boundaries  are  assigned  free- 
slip  conditions,  and  the  outflow  boundary  above  the 
parachute  is  assigned  a  traction- free  condition.  The 
semi-discrete  formulation,  which  is  less  cost-intensive 
than  the  DSD/SST  formulation,  is  adequate  for  the 
stand-alone  simulations  since  there  is  no  time  depen¬ 
dence  in  the  spatial  domain  (i.e.,  no  deformations 
of  the  canopy).  After  this  flow  is  developed,  several 
time  steps  were  computed,  still  with  the  fixed  canopy 
but  by  using  the  DSD/SST  procedure,  to  obtain  the 


starting  CFD  conditions  for  the  FSI  simulation.  The 
developed  pressure  field  and  velocity  field  are  shown 
in  Figure  3. 


Figure  3.  Inflated  T-10:  Initial  CFD  mesh  and  flow- 
field. 

Stand  alone  FSI  simulation : 

Final  conditions  from  the  stand  alone  simula¬ 
tions  are  used  to  initiate  the  FSI  simulation.  Here, 
the  parachute  structure  is  totally  unconstrained.  To 
handle  the  motion  of  the  parachute,  the  fluid  mesh 
moves  globally  relative  to  the  mean  parachute  canopy 
displacement  and  the  canopy  shape  changes  are  ac¬ 
counted  for  in  the  CFD  mesh  using  the  automatic 
mesh  moving  scheme.  The  deforming  parachute, 
with  the  canopy  differential  pressures,  for  the  FSI 
simulation  are  shown  at  four  instants  in  time  in  Fig¬ 
ure  4.  Figure  5  shows  the  computed  drag  force  in 
comparison  to  the  total  gravitational  force  acting  on 
the  parachute  system  (i.e.,  canopy,  suspension  lines, 
risers,  and  payload  weights).  As  expected,  the  drag 
force  oscillates  about  the  weight  of  the  parachute 
system. 

Control  Line  Pulls 

The  ability  to  simulate  riser  pulls  was  initially 
demonstrated  for  a  “soft  landing”  of  a  T-10  FSI  sim¬ 
ulation.13  Here,  a  soft-landing  was  accomplished  by 
smoothly  decreasing  in  time  the  natural  lengths  of 
the  cable  elements  defining  each  of  the  four  risers. 
The  preliminary  soft-landing  simulation  is  being  ex¬ 
tended  to  study  the  response  of  the  parachute  sys¬ 
tem  to  single-  and  double-riser  pulls  and  releases. 
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ulations  for  the  standard  configuration,  for  single- 
and  double-riser  pulls  of  2.5  feet,  and  for  single-  and 
double-riser  releases  of  3.0  feet.  Loss  of  symmetry 
in  the  parachute  canopy  resulting  for  the  riser  pulls 
and  releases  is  evident.  These  asymmetries  provide 
the  canopy  with  lift  which  can  be  used  to  steer  or 
redirect  the  parachute  system. 
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1 — riser  release. 


2 — riser  release. 


Figure  6.  Continued 


tions  for  the  fixed  canopy  configurations.  These  flow 
solutions,  along  with  the  inflated  configurations  for 
the  structure,  will  be  used  as  a  starting  point  for 
future  FSI  simulations. 

The  precise  shape  and  motions  for  the  parachute 
system  with  and  without  control  line  pulls  is  gov¬ 
erned  by  a  strong  interaction  between  the  structure 
and  the  fluid.  Thus,  accurate  prediction  of  riser-slip 
behavior  requires  a  FSI  model.  Results  to  corre¬ 
sponding  FSI  simulations  are  not  yet  available.  Fu¬ 
ture  simulations  will  address  the  FSI  behavior  for 
these  riser-slip  combinations. 

CONCLUDING  REMARKS 


Standard  T-10  (no  contol  line  pulls) 


1-riser  pull  (2.5  feet)  2-riser  pull  (2.5  feet) 


1-riser  release  (3.0  feet)  2-riser  release  (3.0  feet) 

Figure  7.  T-10  canopy  during  standalone  CFD  simu¬ 
lation  (differential  pressures). 


which  can  provide  the  parachute  system  with  limited 
maneuverability.  Riser  slips  and  control  inputs  to 
parachute  systems  (round,  cruciform,  parafoils,  and 
other)  are  real  applications  for  which  this  technol¬ 
ogy  is  being  applied.  The  capability  within  DoD  to 
reduce  the  life-cycle  cost  of  such  systems  should  be 
realized  over  the  next  several  years  as  these  modeling 
capabilities  mature. 


We  have  presented  a  strategy  for  carrying  out 
3-D  simulations  of  parachute  fluid-structure  interac¬ 
tion  has  been  presented  and  demonstrated  for  simu¬ 
lations  of  airdrop  performance  and  control  phenom¬ 
ena  in  terminal  descent.  The  strategy  is  demon¬ 
strated  to  simulate  the  terminal  descent  behavior 
of  a  standard  US  Army  personnel  parachute  system. 
Also,  preliminary  simulations  are  presented  for  the 
response  of  the  parachute  to  a  variety  of  “riser  slips,” 
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